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FOREWORD 

This  report  documents  a  FORTRAN  routine  WOLFQP  for  solving  quadratic 
programming  problems.  Theory,  usage  notes  and  examples  are  included  in  the 
report. 

IRA  M.  BLATSTEIN 
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CHAPTER  1 

WOLFE'S  METHOD  FOR  QUADRATIC  PROGRAMMING 


Wolfe's  quadratic  programming  method^-  is  based  on  solving  a  system  of  linear 
relations  subject  to  complementarity  conditions.  The  linear  system  (without  the 
complementarity  conditions)  can  be  solved  by  the  simplex  method  of  linear 
programming.  The  pivoting  rules  of  the  simplex  method  can  be  restricted  so  that 
the  complementarity  conditions  are  met.  When  the  quadratic  form  in  the  objective 
function  is  positive  definite  (for  a  minimization  problem)  the  pivoting 
restrictions  do  not  prevent  the  simplex  method  from  solving  the  linear  system. 

The  linear  system  is  based  on  the  Karush-Kuhn-Tucker  conditions  (or  Lagrange 
multiplier  rule)  for  the  problem  to  be  solved.  Let  this  problem  be: 


T  T  T1 
minimize  1/2  I X2  J 


Qll  Qi2 

X1 

Ft  tT 

X1 

+[Pl  P2  J 

Ql2T  Q22 

L  -1 

x2 

L—  — 

x2 

+  C 


subject  to  the  constraints 


x^  21  0  (n^  variables)  , 

*2  unconstrained  (n2  variables) , 


[All 

[A21 


<_  b^  (m^  constraints)  , 


|x2 


b£  (ra2  constraints). 


The  constant  term  C  is  included  for  evaluation  purposes  only,  it  does  not 
affect  the  minimizing  point.  If  there  are  no  unconstrained  variables  (n-,  =  0) 
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the  matrix  Q  = 

Q11  Q12 

reduces  to  and  A  = 

fAll  A12] 

reduces  to 

A11 

Q12T  Q22 

<M  ' 
CM 
<1 

pH 

CM 

<  , 

>1_ 

Similar  interpretations  apply  when  n^  =  0,  ni^  =  0  or  =  0.  Inequalities  apply 
componentwise;  i.e.,  x^  >_  0  is  equivalent  to  x^  _>  0,  j  ■  l,...,n^.  Q  is  positive 
definite. 

Introduce  Lagrange  multiplier  vectors  for  the  constraint  x^  0  and  X^ 
and  for  the  other  constraints,  and  form  the  Lagrangian 


L  = 


1/2  [»i 


Q11  Q12 

q12t  q22 


+  c 


+  <  T 
1 


+  ^([hl 
+  X2T([A21 


12 


LX2. 


-b. 


l22l  [*ll  "h) 


from  which,  by  the  Lagrange  multiplier  rule  for  inequality  constrained  problems, 
we  get  the  linear  conditions 


ki+[qii  qi2] 

V 

+  [AuT 

^l1] 

V 

x2 

- 2 

[°12  ^22] 

X1 

+  [V 

A22T] 

*1 

-X2- 

-2- 

+  P1  =  o» 


+  P2  *  0, 


along  with  the  original  constraints 
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X1  i  °» 


the  sign  constraints  on  the  multipliers  for  the  inequalities 


K1  1  °» 

X1  1  °» 

and  the  complementary  slackness  condition 

*bl)+  X2T ([A21  Vfll  -»2)- 

L2. 

Since  each  of  the  n^  +  m^  +  m2  terms  in  this  sum  is  nonpositive  by  the  constraints, 
each  must  individually  vanish: 


Klj  "  °  °r  Xlj  '  °»  j  *  1 . V 

etc. 

An  alternate  phrasing  of  this  system  of  linear  inequalities  and  equations  can 
be  made  if  we  introduce  auxiliary  vector  variables  <2*  and  S2. 


Find  a  solution  to: 


NSWC  TR  82-30 


(1) 


K1  1  0,  <2  =  °»  S1  -  0,  S2  =  0,  x1  >_  0,  \1  >_  0  .  (2) 

and 


(The  complementary  slackness  conditions  (3)  again  require  each  term  in  the 
scalar  products  and  A^  to  vanish.)  If  we  ignore  the  complementary 

slackness  conditions,  we  can  solve  (1)  and  (2)  by  minimizing  a  penalty  function 

ni  n2  mi  m2 

P*  j"!  «ax{0,iclj}+  ^  i  K2j !  +  j=i  maxC-S^  ,0  }  +  ^  |  j 

nl  ®1 

+  ^i  13324  ,0  }  +  max  {-Xjj.0} 

Subject  to  (1)  (and  no  sign  constraints).  The  minimum  value  of  P  is  0  if  and 
only  if  (1)  and  (2)  are  consistent.  Of  course,  it  is  not  necessary  to  weight  all 
the  sign  constraint  violations  equally.  For  example,  we  could  replace  | Sjj |  by 
max{-a  S2 j ,  b  S2j  }  where  a  and  b  are  positive.  Changing  the  weights  in  this 
manner  does  not  affect  the  equivalence  of  the  consistency  of  (1)  and  (2)  to  the 
minimum  being  zero.  A  slightly  different  procedure  is  normally  used  (in  a 
simplex  method  Phase  I).  The  identity  matrix  in  the  first  (n^  +  n£  +  +  n^) 

columns  is  used  as  an  initial  basis  for  the  primal  simplex  method,  any  constraint 
violations  resulting  therefrom  are  penalized,  and  sign  constraints  are  retained 
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in  such  a  way  that  the  penalty  function  is  linear.  For  example,  if  b£j  >  0, 

S2j  *  b2j  in  the  initial  basic  solution.  The  term  •"  +  S2j  +  ***  appears  in  the 
penalty  and  the  restriction  S2j  >  0  is  imposed;  S2j  is  termed  an  "artificial 
variable"  in  this  case.  (If  b2j  were  negative,  the  penalty  would  be  -S2j  and 
S2j  kept  nonpositive.)  For  the  treatment  is  slightly  different.  If  b^j  0, 
set  S-^j  =  b^j,  with  no  penalty.  If  b^j  <  0,  divide  into  the  sum  of  a 
slack  variable  (unpenalized,  nonnegative)  and  an  artificial  variable  S^a 
(penalized,  nonpositive),  set  to  b^j  initially.  This  procedure  is,  however, 
equivalent  to  changing  the  weights  in  P.  Using  b£j  >  0  again  as  an  example,  we 
could  change  the  weighting  so  that  the  penalty  on  S2j  changes  from  | S2j |  to 
max{ (-<*>) S2 j ,  S2j  )  and  in  practice  we  could  replace  ®  by  a  suitable  large  positive 
quantity  M. 

It  is  not  actually  necessary  to  use  a  linear  cost  functional.  With  a  little 
care  in  programming,  the  usual  simplex  method  can  handle  cost  functionals  like  P. 

O 

The  linear  programming  code  LINOPT  uses  the  dual  simplex  method  to  solve  problems 
of  maximizing  a  linear  functional  subject  to  upper  and  lower  bounds  on  all 
variables  and  constrained  quantities.  The  problem  of  minimizing  P  subject  to  (1) 
is  the  dual  of  such  a  problem,  which  could  be  solved  using  LINOPT,  with  the 
effect  of  minimizing  P  subject  to  (1)  by  the  primal  simplex  method. 

Let  us  now  consider  the  complementary  slackness  conditions  (3).  In  terms  of 
the  primal  simplex  method,  they  can  be  enforced  by  ensuring  that  Sjj  and  its 
complementary  variable  X^j  (<^j  an^  xlj)  are  not  simultaneously  basic.  This  is 
an  easily  enforced  pivoting  restriction.  Unfortunately  it  is  too  strong  to  be 
applied  to  minimizing  P  subject  to  (1)  directly:  the  algorithm  may  terminate 
with  P  >  0  even  though  (1),  (2)  and  (3)  are  consistent.  We  get  around  this 
difficulty  by  using  a  phase  I/phase  II  procedure.  Suppose  we  split  into  a 
multiplier  Kj”  ,  subject  to  complementary  pivoting  restrictions,  and  an  artificial 
variable  K^,  not  so  restricted:  +  K*  .  Since  and  cannot  both 

be  basic,  one  of  them  is  always  zero  and  nonbasic. 


We  can  let  be  whichever  of  them  is  basic  (if  one  is)  and  keep  track  of 
whether  this  is  K^a  or  . 
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h 


Suppose  now  that  we  have  a  basic  feasible  solution  to  the  constraints  on  x. 


Then 


I 

0 

A11 

1 

CN 

H 

< 

V 

i— \ 

X) 

1 _ 

0 

I 

A21 

A22_ 

s2 

UJ 

xi 

x2 

S 

1  v 

o,  s 

2  = 

0 ,  x 

1  - 


we  now  minimize  P*  = 


m 


nl 

i  ai  ,  n2 

E 

j=l 

i=i 

1  °. 

xL  >  0,  S1  >_ 

| <2 j I  subject  to  (1)  and 


pivoting  restriction,  the  minimum  will  be  zero,  and  we  will  have  solved  (1),  (2) 
and  (3).  (The  proof  of  this  fact  requires  the  positive  definiteness  of  0.) 


If  we  wish  to  enforce  the  sign  constraints  by  penalty  terms  added  to  P'  these 
terms  must  be  so  heavily  weighted  that  no  change  of  basis  that  reduces  P'  can 
introduce  sign  constraint  violations. 


The  phase  I  is  to  find  the  initial  basic  feasible  solution.  It  is  not 
necessary  to  work  on  a  smaller  problem  using  only  the  [I  A]  part  of  the  full 
coefficient  matrix  fl  M  ^"1  If  we  allow  and  t0  be  basic  and 


[i  o  Q  * 

[o  I  A  0. 


unconstrained  in  phase  I,  they  will  be  basic  throughout  phase  I  and  the  initial 
basis  inverse  for  Phase  II  will  be  the  final  basis  inverse  for  phase  I. 


Phase  I  problem: 


minimize  P, 


M 


mi 

E  max{-S^j  ,0  }  + 
j=l 


m2 

E 

j=l 


nl 


’2j 


+  E  max{-x.  .,0}  +  E  max{-X^.,0 
j=l  2  j=l  2 


subject  to: 
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I 

0 

0 

0 

*11 

Q12 

A  T 
11 

0 

I 

0 

0 

q22 

A  T 
12 

0 

0 

I 

0 

A11 

A12 

0 

0 

0 

0 

I 

A21 

a22 

0 
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V 

A21T 

*2 

= 

-pl 

A22T 

S1 

-p2 

0 

S2 

bl 

0 

X1 

b2_ 

X 

2 

xi 

A  2 

where  M  »  1.  Since  M  scales  out  of  the  problem  it  is  not  actually  needed  here. 
It  is  included  as  part  of  the  setup  for  phase  II.  A  similar  comment  applied  to 
the  penalty  on  A1  which  remains  nonbasic  and  zero  throughout  phase  I  (and  could 
therefore  be  ignored).  Starting  basic  variables  are  <]_,  <2*  S1  and  s2'  Tbf>- 
program  actually  sets  up  and  solves  the  dual  to  this  problem.  Let  k^  be  a  vector 
of  variables  dual  to  k^,  k2  dual  to  k2>^1  dual  to  S^,  dual  to  S,, ,  dual  to 
xl5  dual  to  x2,  dual  to  X1  and  £2  dual  t0  ^2*  For  both  Phases  the  dual 
variables  are  related  by  the  transpose  of  the  coefficient  matrix: 


kl 

*1 

0 

0  0 

k2 

0 

I 

0  0 

ai 

0 

0 

I  0 

ct2 

0 

0 

0  I 

*1 

*12 

A  T  a  F 
11  21 

?2 

Q12T 

Q22 

A  T  A  T 

rtl2  rt22 

*1 

A11 

A12 

0  0 

l2 

A21 

A22 

0  0 

k2 

a2 


(4) 


: 
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The  dual  phase  I  problem  is  to  maximize  z  =  [-p^T  bj_T 

subject  to  the  bounds: 


La2. 


kj  =  0,  k2  =  0,  “M  <  a1  £  0,  -M  <  02  <_  M 
-M  <  q  <  0,  £2  =  0,  -M  <  <  0,  l2  =  0 


(The  bounds  on  and  are  obviously  always  satisfied  in  phase  I.) 

The  objective  function  is  defined  by  the  right-hand  side  of  (1).  The  bounds  are 
defined  by  the  penalties  in  P^:  the  bounds  for  a  dual  variable  deriving  from 
the  penalties  on  the  corresponding  variable.  Violation  of  nonzero  bounds  in 
phases  I  and  II  are  ignored:  they  can  be  made  to  go  away  by  increasing  the  penalty 
on  the  appropriate  variable  and  we  could  have  set  up  the  penalty  function  in  this 
way  in  the  first  place.  Furthermore,  ignoring  violation  of  nonzero  bounds  yields 
a  more  efficient  procedure. 

Assuming  that  phase  I  ends  with  feasibility  shown,  we  start  phase  II  by 
changing  the  penalty  from  P  to 


P  =  P  + 
II  I 


j  =  l 


n? 

nl 

k2J 

j=i 

+  £ 

J-l 

k^  and  k^. 

For  k, 

max  {0,  k”  }, 


-1  ^  k^  £  1,  but  for  k^,  setting  the  bounds  is  complicated  by  the  implicit 
handling  of  the  split  of  into  .  Let  us  start  by  making  each  artificial 

g 

variable  k,  .  basic,  and  treat  k .  as  its  dual,  with  bounds  -1  <  k, .  '  1.  Then  k, . 

lj  lj  -  lj  ~  lj 

will  initially  be  either  +1  or  -1  (depending  on  the  sign  of  ic^).  Ordinarily, 
without  the  pivoting  restriction,  if  k^  =  k^  =  -1,  the  lower  bound  (zero)  on 
k^.  =  k^  =  k^  would  be  violated,  and  we  could  pivot  k^.  into  the  basis  with 
value  zero,  with  no  change  in  the  basis  inverse  matrix.  With  the  pivoting 
restriction,  we  can  still  do  this  if  £  is  not  basic  (and  zero — but  it  will  be 
zero  if  basic  because  we  ignore  violation  of  nonzero  bounds).  Since  we  want  to 
reduce  artificial  variables  to  zero,  we  might  as  well  do  so  at  the  start.  Thus 

the  initial  setup  of  bounds  on  k.  are:  r.  <_  k..  ,  _<  1 

J  j 

where  Tj  =  -1  if  is  basic  at  the  end  of  phase  I, 


Tj  =  0  if  ^  is  nonbasic  at  the  end  of  phase  I. 
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CHAPTER  2 
USER'S  NOTES 


WOLFQP  is  a  FORTRAN  Language  computer  code  written  to  solve  problem  (1)  with 
the  method  described  in  the  previous  section.  The  code  is  a  modified  version 
of  LINOPT  which  takes  advantage  of  the  structure  in  the  Quadratic  Programming 
problem  to  reduce  the  storage  required.  Also  the  user  is  only  required  to  input 
the  QP  problem,  it  is  translated  by  WJLFQP  into  the  Linear  Programming  Problem. ' 

Communication  with  the  calling  program  is  accomplished  through  the  calling 
statement.  The  formal  parameters  are  described  thoroughly  in  the  internal 
documentation  (See  Appendix  A).  This  section  therefore  will  only  attempt  to 
clarify  some  of  the  murkier  details. 

The  scratch  array  SCR  must  be  dimensioned  at  least  N2  +  13N  +  6 
(N  *  N^  +  h’2  +  where  N^  *  Number  of  sign  constrained  variables, 

N2  »  Number  of  sign  unconstrained  variables,  »  Number  of  unequality 
constraints,  M2  =  Number  of  equality  constraints).  Failure  to  do  this  will  cause 
memory  to  be  overwritten  and  result  in  indeterminate  output.  The  scratch  array 
is  divided  into  the  internal  arrays  required  of  LINOPT.  The  correlation  between 
these  internal  variables  and  positions  in  SCR  are  detailed  in  the  internal 
documentation  (See  Appendix  A).  For  detailed  descriptions  of  the  internal 
variables  see  reference  [1]. 


It  Is  noted  that  the  unknown  vector 


X1 

x2 


in  problem  (1)  has  sign  constrained 


xj  and  sign  unconstrained  X2  variables.  For  ease  of  notation  the  constrained  and 
free  variables  are  grouped  separately.  WOLFQP  does  not  require  the  user  to  so 
order  the  variables.  The  integer  vector  IND  allows  the  user  to  permute  the 
variables  without  permuting  the  Q,  A1  and  A2  matrices.  This  allows  the  user 
to  form  the  matrices  in  a  convenient  manner. 

The  First  N^  locations  of  IND  contain  the  indices  for  which  there  are  sign 
constraints.  Locations  +  1...,  N]_  +  N2  contain  the  indices  for  the  free 
variables. 
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Also  it  may  be  inconvenient  to  reform  the  constraint  matrices  for  a  call  to 
WOLFQP.  As  an  example  the  user  may  wish  to  call  WOLFQP  with  only  some  of  the 
constraint  rows  active.  Locations  Np  +  N2  +  l....N]_  +  N2  +  M^  of  IND  contain 
the  indices  of  the  active  rows,  if  any,  in  the  inequality  constraint  matrix  Al. 

In  the  same  manner  locations  Np  +  N2  +  +  I...N1  +  N2  +  Ml  +  M2  of  IND  contain 

the  indices  of  the  active  rows,  if  any,  in  A2;  the  equality  constraint  matrix. 

If  the  value  of  the  objective  function  is  desired  the  user  should  set  the 
formal  parameter  VALUE  to  a  non-zero  value.  This  results  in  the  calculation 
of  the  objectives  function  value  at  the  optimal  solution.  The  Scalar  C  is  ignored 
except  in  this  calculation. 

Two  inputs  which  control  the  algorithm  are  ITMAX  and  NINVT,  corresponding 
to  IPASS  (2)  and  IPASS(8)  respectively.  ITMAX  is  the  iteration  limit.  Upon 
completion  of  ITMAX  iterations  control  is  returned  to  the  calling  program  with 
IERR  (IPASS  (4))  set  to  2.  All  internal  storage  has  been  saved.  If  desired 
WOLFQP  may  be  recalled  and  computation  continued  by  calling  WOLFQP  with  FLAG 
set  to  true.  No  other  change  is  necessary  or  advised. 

On  output  the  optimal  distribution  is  placed  in  the  array  Y.  The  value 
is  calculated,  if  desired.  The  iteration  number  and  error  indicator  are  stored 
in  the  appropriate  positions  of  IPASS  (See  Appendix  A).  Also  if  the  user  desires 
to  look  at  the  internal  variables  SCR  contains  these  values. 

Other  error  conditions  (IERR  =  1,3,4+)  must  be  corrected  before  re-calling 
WOLFQP.  For  these  conditions  the  tableau  must  be  re- initialized.  (Flag  set  to 
False) . 

An  upper  limit  on  ITMAX  should  be  between  5*N  and  10*N.  Some  problems  may 
take  more;  most  should  take  fewer  iterations.  Practically  ITMAX  should  be  set  so 
that  an  inordinate  amount  of  CP  time  is  not  consumed  before  the  problem  formulation 
has  been  thoroughly  checked. 

NINVT  controls  the  re-inversion  of  the  tableau.  Since  succeeding  inverses 
are  formed  as  perturbations  of  previous  inverses  truncation  error  can 
accumulate  in  the  inverse.  To  remedy  this  situation  the  inverse  must  be 
re-formed.  Every  NINVT  iterations  the  inverse  will  be  reformed.  To  avoid 
unnecessary  calculations  this  should  not  occur  too  frequently.  The  recommended 
value  is  between  2  x  N  and  3  x  N. 

The  final  input  controlling  calculations  is  EPSl.  This  is  used  to  control 
zero  tolerance.  Zero  tolerance  is  important  in  two  parts  of  the  algorithm. 

Tableau  entries  which  are  less  than  EPSl  in  magnitude  are  treated  as  zero. 
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Their  values  are  considered  to  be  truncation  error.  Thus  for  this  purpose  EPS1 
is  set  to  the  truncation  error  of  the  machine. 

The  other  spot  where  zero  tolerance  is  involved  is  when  constraint  violations 
are  checked.  Constraint  violations  less  than  EPS1  are  treated  as  zero.  This 
is  also  to  avoid  truncation  error  but  the  user  must  be  careful.  If  the  problem 
is  scaled  such  that  entries  in  the  tableau  are  less  than  EPS1  their  entries 
will  be  ignored.  For  this  purpose  EPS1  must  be  on  the  order  of  the  smallest 
entries  in  the  tableau  (Q,  Al  and  A2).  The  input  value  for  EPS1  should  be  the 
minimum  of  these  two  values. 

The  major  modification  to  LINOPT  occurs  in  PIVROW.  This  subroutine  chooses 
the  incoming  basic  variable  which  is  required  to  meet  the  complementary  slackness 
condition.  Flags  contained  in  the  internal  array  BASIC  are  used  to  facilitate 
this  checking. 

Each  word  in  the  BASIC  array  contains  seven  bits  used  as  flags  and  if  the 
complementary  slackness  condition  applied  to  the  variable,  (a  complementary  index). 
Bits  are  used  to  reduce  the  storage  requirement.  The  layout  of  each  basic  word 
is  described  in  the  internal  documentation  of  SETQP.  (See  Appendix  A). 

The  Bit  operations  OR,  AND,  SHIFT  and  COMPLEMENT  are  required  to  access 
and  define  the  flag  bits.  These  operations  are  described  more  thoroughly  in 
machine  dependencies  (Appendix  C) . 

The  Flag  Bits  are  also  used  in  DSIMP.  In  DS IMP  the  basic  or  non-basic  bit  is 
set  as  appropriate  and  bounds  are  reset,  if  applicable,  on  variables  leaving 
the  basis. 
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CHAPTER  3 
EXAMPLES 

A  short  ;est  program  appears  at  the  end  of  the  source  listing.  (Appendix  A). 
The  problem  is  defined  by  the  NAMELIST  inputs.  Mnemonic  names  are  used  in  the 
NAMELIST  which  are  then  loaded  into  the  IPASS  array.  The  following  examples 
list  the  NAMELIST  inputs  and  the  output  resulting  from  a  subsequent  call  to  QPTAB. 
Example  1. 

This  problem  illustrates  the  basic  use  of  the  algorithm. 


min  1/2  [xj_,  x2] 

4  -2 

X1 

+  [6  0] 

X1 

-2  4 

—  — 

x2 

x2 

xi,  x2  >  0 
xi  <  2 

x2  <  1 
xL  +  x2  =  2 

This  problem  requires  the  following  inputs 
$IN 

EPS  =  l.E-10, 

FLAG  =  .FALSE., 

IND  =  1,2, 1,2,1, 

N  =  5, 

NCON  =  2, 

NUNC  =  0, 

NINEQ  =  2, 

NINV  *  3000, 

ITMAX  =  100, 

Q  =  4. ,  -2. , 

Q (1 , 2)  =  2. ,4., 

A1  =>  1.,  0., 

Al(l, 2)  =  0. ,1. , 

A2(l,l)  =  1., 

A2(l,2)  =  1., 

OBJ  =  6. ,0. ,2. ,1. ,2. , 

C  -  0.  , 

VALUE  =1., 

SEND 
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output  of  QPTAB  after  call  to  WOLFQP 
$0UT 

N  =5, 

NCON  =  2 , 

NUNC  =  0, 

NINEQ  =  2, 

NEQ  =  1, 

NINV  =  3000, 

ITMAX  =  100, 

ITER  =  2, 

IERR  =  0, 

C  =  0.0, 

VALUE  =  . 8E+01, 

$END 

Q 

. 4000E+01  2000E+01 

2000E+01  . 4000E+01 

P 

. 6000E+01  0. 

A1*X  B1  A1 

. 1000E+01  . 2000E+01  .1000E+01  0. 

. 1000E+01  . 1000E+01  0.  . 1000E+01 

A2*X  B2  A2 

. 2000E+01  . 2000E+01  .1000E+01  .1000E+01 

X 

. 1000E+01  . 1000E+01 

Basically  the  output  reprints  the  inputs  with  the  addition  of  the  solution 
vector  X  =  [l.,l.],  the  value  of  the  function  VALUE  =  8.,  (VALUE  on  input  was 
non-zero).  The  number  of  iterations  required  ITER  =  2,  the  error  indicator 
IERR = 0  indicating  solution  and  each  constraint  is  evaluated  as  a  check. 

(See  A1*X  and  A2*X  columns.  • 
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Example  2. 

This  example  illustrates  how  sign  constraints  can  be  written  implicitly 
or  explicitly.  The  inputs  differ  slightly.  The  explicit  version  requires  more 

iterations. 

min  1/2  [xj^,  x2,  X3,  X4] 


■l  0  6  1  " 

-  i 

X1 

0  4  6  4 

x2 

6  6  61  12 

x3 

1  4  12  100 

*4 

«  - 

+  [0, 0, 0,  - 1] 


‘X1 

x2 

x- 


SL’  x2> 


x4  1  0 


Inputs  for  Problem  1 
$IN 

EPS  =  l.E-10, 

FLAG  =  .FALSE. , 

OBJ  =  10*0. , 

N  *  4.  , 

NCON  =  4, 

NUNC  *  0, 

NINEQ  =  0, 

NINV  -  2000, 

ITMAX  =  100, 

0(1,1)  -  1. ,0. ,6. ,1. , 

Q(l»2)  =  0. ,4. ,6. ,4. , 

Q(l, 3)  “  6 • , 6 • , 61 • , 12 • , 

Q( 1,4)  *  1 . , 4 . , 12 . , 100 . , 

OBJ  =  0.,0.,0.,-l., 

C  *  0.  , 

VALUE  =  1. , 

IND  =  1,2, 3, 4, 

$END 

One  notes  that  the  sign  constraints  are  entered  implicitly  and  are  the  only 
constraints. 
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Outputs : 


$OUT 

N  »  4, 

NCON  -  4 , 

NUNC  -  0, 

NINEQ  -  0, 

NEQ  =  0, 

NINV  =  2000 , 

ITMAX  =  100, 

ITER  -  1, 


IERR  =  0, 

C  =  0.0, 


VALUE 

II 

1 

(. 

NJ 

\  i 

$END 

Q 

;  1 

1 

. 1000E+01  0. 

0.  . 4000E+01 

. 6000E+01  • 6000E+01 

. 1000E+01  . 4000E+01 

. 6000E+01 
. 6000E+01 
. 6100E+02 
. 1200E+02 

. 1000E+01 
. 4000E+01 
. 1200E+02 
. 1000E+03 

! 

;  1 

0. 

P 

0. 

0. 

-. 1000E+01 

l 

0. 

X 

0. 

0. 

. 1000E-01 

Note  only  1  Iteration  is  required. 

The  following  inputs  solve  the  same  problem,  but  x-^  and  x2  are  implicitly 
constrained.  X3  and  X4  have  explicit  constraints  in  the  A1  matrix. 

Inputs : 

$IN 

EPS  =  l.E-10, 

FLAG  =  .FALSE., 

A1  -  300*0.,  A2  =  300*0.,  Q=900*0.,  OBJ  =  30*0., 

OBJ  =  10*0, 

N  =  6, 

NCON  =2, 

NUNC  =  2, 

NINEQ  =  2, 

NINV  *=  2000, 

ITMAX  -  100,  * 

Q(l,l)  “  1. ,0. ,6. ,1. , 
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Q(l,2)  =  0.,4.,6.,4., 
Q(l,3)  a  6. ,6. ,61. ,12. , 
Q(l,4)  ■  1. ,4. ,12. ,100. , 
Al(l,l)  =  -1 • , 0. ,0 . , 0. , 
Al(l,2)  *  0 . 1. ,0. ,0 . , 
Al(l,3)  =  0. ,0. 1. ,0. , 
Al(l,4)  =  0. ,0. ,0. ,— 1. , 
OBJ*  *  0 . , 0 • , 0. ,— 1 . , 

C  =  0. , 

VALUE  -  1., 

IND  =  1,2, 3, 4, 3, 4, 

$END 


Output : 


$OUT 

N 

= 

6, 

NCON 

= 

2, 

NUNC 

= 

2, 

NINEQ 

= 

2, 

NEQ 

= 

o. 

NINV 

= 

2000, 

ITMAX 

= 

100, 

ITER 

s 

3, 

IERR 

- 

o, 

C 

* 

0.0, 

VALUE 

$END 

s 

5E-02 

Q 

. 1000E+01 

0. 

. 6000E+01 
. 1000E+01 


0. 

.4000E+01 
. 6000E+01 
.4000E+01 


.  6000E+01 
.  6000E+01 
.6100E+02 
. 1200E+02 


. 1000E+01 
. 4000E+01 
. 1200E+02 
. 1000E+03 


P 

0.  0. 

A1*X  B1 

0.  0. 

- . 1000E-01  0. 


0.  1000E+01 

A1 

0.  0.  1000E+01  0. 

q  0.  0.  -.1000E+01 


0. 


X 


0. 


0. 


. 1000E-01 


Note  the  solution  is  the  same  as  the  implicitly  constrained  example  but  this 
solution  required  3  interations  instead  of  1  the  problem  has  increased  in  size 
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from  N  =  4  to  N  =  6  thus  requiring  additional  storage. 

If  one  wishes  to  continue  in  this  vein,  all  sign  constraints  can  be  made 
explicit.  The  dimension  of  the  problem  becomes  N  “  8.  The  same  solution  now 
requires  7  iterations.  The  reader  may  verify  this  if  desired. 

Example  3 

This  is  an  example  displaying  the  restart  capabilities  of  the  code.  Since 
the  Q  matrix  is  the  identity  this  is  a  constrained  least  squares  problem  of  a 
particularly  simple  type.  By  following  the  directions  in  Appendix  B  the  user  can 
eliminate  the  need  to  store  the  identity  matrix  in  memory. 

This  problem  arises  as  a  sub-problem  of  the  feasible  direction  method.  What 
is  required  is  to  find  the  direction  of  steepest  descent  subject  to  the  binding  ' 
constraints.  The  feasible  direction  method  can  at  times  require  many  constrained 
gradients  to  be  calculated.  The  constraint  matrix  is  usually  full  throughout 
the  process  but  which  constraints  are  binding  depend  on  the  given  position  where 
the  constrained  gradient  is  to  be  calculated. 

This  motivates  the  use  of  the  linked  list  IND.  First  to  indicate  in  the  first 
N^  positions  the  indices  of  the  sign  constrained  variables.  In  the  following 
N2  positions  the  indices  of  sign  unconstrained  variables.  Thus  the  variables 
need  not  be  ordered  constrained  followed  by  unconstrained.  The  next  positions 
of  IND  are  filled  with  the  row  indices  of  the  active  inequality  constraints  and 
finally  the  last  M2  positions  are  filled  by  the  indices  of  the  active  equality 
constraints. 

Thus  for  the  feasible  direction  method  the  full  constraint  matrix  is  formed 
once.  Afterward  only  the  row  indices  of  the  matrices  must  be  manipulated  to 
pass  the  current  constraint  set.  Note  the  objective  row  containing  the  right 
hand  side  of  the  constraints  must  be  reformed  each  time  it  changes. 

The  gradient  projection  problem  can  be  stated: 

min  1/2  <  g  -  g(x) ,  g  -  g(x)  > 
g 

where  g(x)  is  the  gradient  at  x  and  g  is  the  constrained  gradient. 

Subject  to 

0  for  all  x^  =  0 
Maintaining  x^  as  non-negative 
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A1  g  1  0 
A2  g  -  o 

where  for  the  gradient  project  problem  and  A2  are  the  constraints  matrices  on 

K. 


Inputs 

$IN 

EPS  *  l.E-10 
FLAG  *  .FALSE., 

C  -  0.  , 

VALUE  *  1. , 

IT  -  26, 

NCON  -  12, 

NUNC  =  8, 

NINEQ  «*  4, 

NINV  =*  3000, 

ITMAX  -  100, 

Q  =  900  *  0.  , 

Q(l»l)  -  1., 

Q(2,2)  =  1., 

Q(3,3)  =  1., 

Q(4,4)  =  1., 

Q(5,5)  -  1., 

Q(6,6)  =  1., 

Q(7,7)  =  1., 

Q(8,8)  *  1., 

Q(9,9)  =  1., 

Q(10, 10)  -  1., 

Q(ll.ll)  -  1.. 

Q(12, 12)  *  1., 
Q(13,13)  *  1., 
Q(14,14)  =  1., 
Q(15,15)  *  1., 

Q(16, 16)  *  1., 

Q(17 , 17)  -  1., 
Q(18,18)  =  1., 
Q(19,19)  »  1., 
Q(20,20)  =  1., 

A1  =  300*0., 

Al(l,l)  =  1., 

Al(l, 2)  -  1., 

Al(l, 3)  =  1., 

Al(l, 4)  =  0.  ,1. , 
Al(l, 5)  -  0. , 1.  , 
Al(l, 6)  =  0.,1., 
Al(l, 11)  =  2*0., 1., 
Al(l, 12)  *  2*0. ,1., 
Al(l, 13)  *  2*0., 1., 
Al(l, 14)  -  3*0. ,1., 
Al(l, 15)  =  3*0. ,1., 
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Al(l, 16)  =  3*0., 1., 

A2 (1, 1)  =  300*0.  , 

A2(l,l)  *  1., 

A2(l,2)  -  1.  , 

A2(l,3)  =  1., 

A2 (1,4)  =  1., 

A2(l ,  5)  =  1., 

A2(l,6)  =  1., 

A2 ( 1 , 7 )  =  1., 

A2 (1, 8)  =  1., 

A2 (1, 9)  =  1., 

A2(l,10)  =  1., 

A2 (1,11)  =  0.,1., 

A2(l,12)  =  0. ,1.  , 

A2(l,13)  =  0.,1., 

A2(l,14)  =  0.,1., 

A2 (1,15)  =  0.,1., 

A2(l,16)  =  0. ,1. , 

A2(l,17)  =  0. ,1. , 

A2(l,18)  =  0.,1., 

A2(l,19)  =  0.,1., 

A2 (1,20)  =  0.,1., 

OBJ  =  -1. ,-.566, -.8187, -.9927, -.8187, -.8187, -.8047, -.7800, -.9643, 
-.7701, -.0659, -.0221, -.0130, -.0667, -.0158, -.0043, 14*0. , 

IND  =  2,3,5,6,8,10,11,12,13,18,19,20,1,4,7,9,14,15,16,17,1,2,3,4,1,2, 
$END 


Outputs: 

tour 


Outputs  Continued: 
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Note  the  constraint  matrices.  The  unequality  constraints  are  groups  of 
three  consecutive  variables.  The  equality  constraints  are  groups  of  ten 
consecutive  variables.  Also  26  iterations  were  required  to  determine  the 
solution. 

This  next  problem  is  identical  to  the  previous  but  instead  of  using  the 
permutation  indice  the  matrices  have  been  reformed.  Sign  constrained  variables 
have  been  ordered  first  followed  by  sign  unconstrained  variables. 

If  this  must  be  done  many  times  for  a  given  constraint  tableau  one  can  see 
a  considerable  work  would  be  required.  Also  ITMAX  has  been  set  to  5. 


Inputs: 

$IN 

EPS  =  l.E-10, 

FLAG  =  .FALSE. , 

IND  -  1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,1,2,3,4,1,2, 
C  =  0., 

VALUE  «  1., 

N  *  26, 

NCON  =  12, 

NUNC  =  8, 

NINEQ  *  4, 

NINV  =  3000, 

ITMAX  *  5, 

Q  *  900  *  0. , 

Q(l,l)  =  1., 

Q(  2,2)  =  1., 

Q(3,3)  -  1., 

Q(4,4)  =  1., 

Q<5,5)  =  1., 

Q(6,6)  =  1., 

Q(7,7)  =  1., 

Q(8,8)  =  1., 

Q(9,9)  *  1., 

Q(10, 10)  =  1., 

0(11,11)  =  1., 

0(12,12)  =  1., 

0(13,13)  =  1., 

0(14,14)  =  1., 

0(15,15)  =  1., 

0(16,16)  =  1., 

0(17,17)  =  1., 

0(18,18)  *  1., 

0(19,19)  =  1., 

0(20,20)  *  1., 

A1  =  300*0. , 

Al(l,l)  =  1., 

Al(l, 2)  =  1., 

Al(l, 3)  =  0.,1., 
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Al(l, 4)  -  0.  ,1. , 

Al(l, 7)  =  2*0. ,1., 

Al(l, 8)  «  2*0., 1., 

Al(l, 9)  =  2*0., 1., 

Al(l, 13)  =  1., 

Al(l, 14)  =  0..1., 

Al(l, 17)  =  3*0., 1., 

Al(l, 18)  -  3*0. ,1., 

Al(l, 19)  =  3*0., 0., 

A2(l,l)  =  300*0., 

A2(l,l)  =  1., 

A2(l,2)  -  1., 

A2 ( 1 , 3 )  =  1., 

A2(l,4)  =  1., 

A2 (1, 5)  =  1., 

A2(l,6)  =  1., 

A2 (1, 7)  *  0. ,1. , 

A2(l,8)  =  0.,1., 

A2  (1,9)  =  0.,1., 

A2 (1,10)  =  0.,1., 

A2(l,ll)  =  0.,1., 

A2 (1, 12)  =  0.,1., 

A2 (1,13)  =  1., 

A2(l,14)  =  1., 

A2  (1,15)  =  1., 

A2 (1 , 16)  =  1., 

A2 (1,17)  *  0..1., 

A2(l,18)  -  0.,1., 

A2(l,19)  =  0..1., 

A2(l, 20)  =  0.,1., 

OBJ  =  -.8566, -.8187, -.8187, -.8187, -.78, -.7701, -.0659 ,-.0221, -.0130, 
3*0., -1. ,-.9927, -.8047, -.9643, -.0667, -.0158, -.0043, 11*0. , 

$END 
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Outputs 

SOUl 


1030  0  3  0  0  9  1 


ooaooooeoooooo 


.0330000000®  ooo: 


»OOOOOOOOOOOC®  o  o  o  o  I 


13000000  30303< 


100000030  0003000 


o  o  o  o  o  o  o 


H  II  N  H  «  tl  N  II  n  n  N  , 


9  *< 

Z  u  w  >  <  a  a 

ozzcizziiz 

<JD*"t*J*"»*»*“*u 


13  0  0  0  0  3 


Output  Continued: 
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Note  IERR  returns  equal  to  2  indicating  the  iteration  limit  has  been  reached. 
To  continue  we  call  WOLFQP  changing  MAXIT  to  100.  No  other  change  is  necessary. 


Inputs: 


$IN 

EPS  *  l.E-10, 

FLAG  =  .TRUE., 

IND  »  1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,1,2,3,4,1,2, 
C  =  0., 

VALUE  =  1., 

N  =  26, 

NCON  =  12, 

NUNC  =  8, 

NINEQ  =  4, 

NINV  =  3000, 

ITMAX  =  100, 

Q  *  900  *  0. , 

Q(l,l)  =  1., 

Q(2,2)  *  1., 

Q(3, 3)  -  1., 

Q(4,4)  =  1., 

0(5,5)  -  1., 

0(6,6)  -  1., 

0(7,7)  =  1., 

0(8,8)  *  1., 

0(9,9)  =  1., 

0(10,10)  =  1., 

0(11,11)  =  1., 

0(12,12)  =  1., 

Q(13,13)  *  1., 

0(14,14)  =  1., 

0(15,15)  *  1., 

0(16,16)  =  1., 

Q(17,17)  =  1., 

Q(18,18)  =  1., 

Q(19,19)  =  1., 

Q(20,20)  =  1. , 

A1  -  300*0., 

Al(l.l)  =  1.  , 

Al(l, 2)  =  1., 

Al(l, 3)  =  0. ,1. , 

Al(l, 4)  =  0.,1., 

Al(l, 7)  =  2*0., 1., 

Al(l, 8)  =  2*0. ,1., 

Al(l, 9)  =  2*0. ,1., 

Al(l, 13)  =  1., 

Al(l, 14)  =  0.,1., 

Al(l, 17)  *  3*0., 1., 

Al(l, 18)  =  3*0. ,1. , 

Al(l, 19)  -  3*0. ,1., 

A2  (1,1)  =  300*0., 

A2 (1,1)  =  1., 
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A2(l,2)  »  1., 

A2 (1 , 3)  -  1., 

A2 (1,4)  -  1., 

A2(l,5)  -  1., 

A2(l,6)  -  1., 

A2 ( 1 , 7 )  -  0.,1., 

A2 ( 1 , 8 )  *  0.,1. , 

A2(l,9)  -  0. ,1. , 

A2 (1,10)  -  0.,1., 

A2(l,ll)  -  0..1., 

A2(l,12)  -  0.,1., 

A2(l,13)  -  1., 

A2(l,14)  -  1., 

A2(l,15)  -  1., 

A2(l,16)  -  1., 

A2(l,17)  -  0.,1., 

A2(l,18)  -  0.,1., 

A2(l,19)  *  0..1-, 

A2(l,20)  =  0.,1., 

OBJ  =  -.8566, -.8187, -.8187, -8187, -.78, -.7701, -.0659, -.0221, -.0130, 
3*0. ,-l. , -.9927, -.8047, -.9643, -.0667, -.0158, -.0043, 11*0. , 

END 


« 
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After  unscrambling  the  output  so  that  it  corresponds  to  the  first  problem 
of  this  example  set  one  can  see  the  answers  are  the  same. 

Example  4 

Example  3  illustrated  a  constrained  least-squares  problem  which  consisted 
of  finding  a  Euclidean  projection  onto  a  cone.  More  general  least-squares 
approximation  problems  lead,  by  the  same  algebraic  manipulations  which  yield 
the  normal  equations,  to  quadratic  minimization  problems  with  the  Q  -  matrix 
not  equal  to  the  identity. 

Given  k  values  Hi,...,  H^,  we  wish  to  find  n  values  x^ .  Xjj  which  minimize 

k  2  ri 

1/2  h.  where  h^  *  Xj  -  H.,  the  i-th  residual,  subject  to  upper 

i=l  1  i=l  ■LJ  J  i 

and  lower  bounds  on  the  variables  x j ,  j»l,...,n.  In  matrix  terms  we  can  define 
a  vector  residual  h  *  CX  -  H,  where  C  iskx  n.  Expanding  the  quadratic  gives 
us  the  problem: 

minimize  1/2  xT(CTC)X  -(HTC)x  +  1/2  HTH 


subject  to  [  X]  X  <  [U  j 

T 

where  u  and  Z  are  vectors  of  upper  and  lower  bounds.  Here  C  C  corresponds  to 
Q,  -HtC  to  p  and  1/2  HTH  to  C.  (More  general  linear  constraints  could  also 
be  used.) 

The  following  example,  for  which  we  are  indebted  to  George  Gray  of  E22, 

NSWC,  treats  the  approximation  of  a  given  magnetic  field.  The  original  data, 
unsealed,  make  all  the  Q-matrix  entries  very  small  (~  10"^).  Since  values  less 
than  EPS1  in  magnitude  are  set  to  zero  in  the  tableau  calculations,  we  must 
ensure  that  EPS1  is  small  enough  to  ensure  that  no  significant  Q  -  matrix  entries 
are  zeroed.  We  present  three  cases,  one  with  EPS1  =  10"^-®  (which  bombs),  one 
with  EPS1  *  0,  and  one  with  EPS1  =  10-^®  but  with  the  data  rescaled  (x  divided 
by  10^,  C  multiplied  by  10^)  so  that  this  difficulty  does  not  arise. 
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APPENDIX  B 

MODIFICATION  FOR  PROJECTION  PROBLEMS 


WOLFQP  as  written  will  solve  a  projection  problem.  For  large  problems  space 
may  be  reduced  by  modification  to  the  program. 

In  the  projection  problem  the  Q  matrix  reduces  to  the  identity.  This  need 
not  be  stored.  The  specific  modifications  are: 

WOLFQP 

Delete  the  Do  Loop  with  130 
Replace  with 

TEMP  =  TEMP  +  SCR(IP3  +  I  +  N) 

GETROW 

Delete  Do  Loop  with  50 
Replace  with 

AROW  (J)  =  AROW(J)  +  E(KK  +  JJ) 

PSOL 

Delete  Do  Loop  with  40 
Replace  with 

X(KI)  =  X(KI)  +  X(KK) 

By  then  substituting  a  dummy  argument  for  Q  in  the  calling  Sequence  the  user 
has  removed  the  necessity  of  storing  the  identity  matrix.  Of  course  the  user 
may  then  remove  Q  and  NQD  completely  from  the  calling  sequence  if  desired. 
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APPENDIX  C 

MACHINE  DEPENDENCIES 


BIGM  -  Defined  by  Data  Statement  in  WOLFQP.  BIGM  represents 

machine  infinity.  For  CDC  Machines  10^0  is  used. 

This  constant  clearly  depends  on  the  exponent  range 
available. 

BASIC  -  An  array  used  to  store  internal  flags.  To  reduce  space 

7  bits  of  each  BASIC  word  are  used  to  store  flags  the 
remaining  bits  are  required  to  store  an  integer.  If 
the  word  size  of  the  machine  is  M  bits  this  integer  must 
be  less  than  2^~  -1.  For  CDC  machines  this  becomes  2^^-i. 
This  number  2M  -1  represents  the  largest  problem  which 
can  be  solved  with  this  encoding  of  the  algorithm. 

Since  bit  operations  are  performed  on  the  Basic  array  certain  bit  functions 
must  be  utilized. 

Octal  constants  are  used  to  set  certain  bit  patterns  in  SETQP,  DSIMP  and 
PIVROW.  They  are  either  used  to  set  the  bit  patterns  or  mask  certain  bits  from 
the  given  word.  The  octal  constants  appear  only  in  Data  Statements  in  the 
appropriate  routines.  They  appear  as  integers  followed  by  the  letter  B.  These 
are  numbers  written  in  Base  8  and  must  be  written  with  the  appropriate 
notation  and  in  the  appropriate  base  for  the  machine. 

i.e.  7B  =  7  Hex 

10B  =  8  Hex 

60B  =  30  Hex 

The  Machine  functions  utilized  are: 

0R(A,B)  -  Perform  bit  by  bit  logical  or 

on  A,B.  The  truth  table  is 

0 _ 1_ 

0  0  1 

111 


95 


N'SWC  TR  82-30 


AND(A,B)  -  Perform  bit  by  bit  logical  and 
on  A,B.  The  truth  table  is 


,_0 _ 1_ 

0  0  0 
10  1 

SHIFT(A,I)  -  SHIFT  the  word  I  bits.  If  I  is  positive  SHIFT  left  I  positions. 
If  I  is  negative  shift  right  I  positions, 
i.e.  if  A  contains  the  following  bit  pattern 

A  =  0001011010 

then 

SHIFT(A,2)  =  0101101000 
SHIFT (A, -2 )=  0000010110 


COMPL(A)  -  form  bit  by  bit  Boolean  complement  of  A. 
i.e.  if  A  =  10110 
COMPL(A)=  01001 

The  following  library  function  is  required  in  PIVCOL 

RANF(X)  -  returns  a  uniform  random  number,  the  argument  is  ignored,  in 
the  function. 

The  random  number  is  used  in  an  anti-cycling  procedure.  Some 
other  method  of  preventing  cycling  could  be  used  instead. 
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